SCIENTIFIC 

REPORTS 



s 




SUBJECT AREAS: 

MACROECOLOGY 

ASTEROIDS, COMETS AND 
KUIPER BELT 

PALAEONTOLOGY 

ECOLOGICAL NETWORKS 

Received 
6 February 201 3 



Accepted 
8 April 201 3 

Published 
7 May 2013 



Correspondence and 
requests for materials 
should be addressed to 
DA.V. (daril@uw.edu) 
or G.P.W. (gpwilson@ 
u. washington.edu) 



Bivalve netv\/ork reveals latitudinal 
selectivity gradient at the end-Cretaceous 
mass extinction 

Daril A. Vilhena\ Elisha B. Harris^ Carl T. Bergstrom^ ^ Max E. Maliska\ Peter D. Ward^-^ ^ 
Christian A. Sidor^-^, Caroline A. E. Stromberg^-^ & Gregory P. Wilson^ 

^ Department of Biology, University of Washington, Seattle, WA 98 1 95-1 800, ■^Burke Museum of Natural History and Culture, 
University of Washington, Seattle, WA 981 95-301 0, ^Santa Fe Institute, 1 399 Hyde Park Rd., Santa Fe, NM 87501 , "^Department 
of Earth and Space Science, University of Washington, Seattle, WA 98 1 95-1 800. 

Biogeographic patterns of survival help constrain the causal factors responsible for mass extinction. To test 
whether biogeography influenced end-Cretaceous (K-Pg) extinction patterns, we used a network approach 
to delimit biogeographic units (BUs) above the species level in a global Maastrichtian database of 329 bivalve 
genera. Geographic range is thought to buffer taxa from extinction, but the number of BUs a taxon occurred 
in superseded geographic range as an extinction predictor. Geographically, we found a latitudinal selectivity 
gradient for geographic range in the K-Pg, such that higher latitude BUs had lower extinction than expected 
given the geographic ranges of the genera, implying that (i) high latitude BUs were more resistant to 
extinction, (ii) the intensity of the K-Pg kill mechanism declined with distance from the tropics, or (iii) both. 
Our results highlight the importance of macroecological structure in constraining causal mechanisms of 
extinction and estimating extinction risk of taxa. 

Mass extinctions have disproportionately shaped the evolutionary history of life\ During these rare, 
geologically rapid events, the rules of selectivity that prevail in background extinctions do not always 
clearly apply^. To delimit what may be independent selection processes and constrain possible causal 
factors, paleontologists have sought biologically meaningful patterns of survivorship in mass extinctions^'^. In the 
marine realm, there are some taxon -specific cases where survivorship is linked to ecological traits - for example, 
in the K-Pg mass extinction event, reliance on photosymbiosis among scleractinian corals severely reduced 
survivorships and sea urchin feeding strategy correlates positively with survivorships. However, more often, 
survivorship in mass extinctions appears to be linked to increased geographic range size above the species leveF, 
suggesting that biogeographic history of a taxon plays a vital role in sorting survivors from victims. 

The correlation of geographic range at the lineage level with lineage survival, but apparent lack of physiological 
mechanisms in determining survivorship, suggests that lineages with similar biogeographic histories should have 
similar chances to survive mass extinction. Histories of ocean surface currents, plate tectonics, and environments 
shape biogeographic patterns above the species level by allowing lineages (genera) dispersal opportunities. Here, 
we explore whether these emergent biogeographic regions impact our understanding of extinction processes. 

Here we employ a network approach to test whether biogeographic structure above the species leveP"^^ 
correlates with bivalve survivorship in the K-Pg. Network methods have been useful in a broad range of applica- 
tions, for example, to model the transmission of disease in social networks to describe the structure of scholarly 
communication^^, and to model the stability of ecosystems in response to extinction^^. A network approach can 
also reveal spatial patterns of taxa from geographic range data^^. We generated a Maastrichtian network from the 
bivalve dataset, and used stratigraphic ranges from the Paleobiology Database to determine which genera 
survived the K-Pg mass extinction event^^. We adapt a network-based clustering approach^^ to reveal biogeo- 
graphic units (BUs) from patterns of geographic ranges (see methods). 

The K-Pg event (ca. 66 Ma) is an ideal case to test whether biogeography influences survival in mass extinc- 
tions. It was geologically abrupt and was associated with significant changes to marine productivity and ocean 
chemistry, dramatic restructuring of marine and terrestrial communities, long-term effects on evolutionary rates 
and biogeography^ and the extinction of up to 76% of all species^^. As the most recent of the "big five" mass 
extinctions, the quantity, quality, and spatial resolution of the geological and paleontological data for the K-Pg 
interval are also better than those available for more ancient mass extinctions^^. Bivalves have emerged as a model 
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system for examining macroevolutionary phenomena due to their 
excellent preservational record^^, deep evolutionary history, spatial 
ubiquity, and the significant effort dedicated to standardizing their 
systematics^. 

The bivalve dataset that was used in this study was downloaded 
from the Paleobiology Database (PBDB^^,) and consists of 3,445 
occurrences of 329 bivalve genera from 105 Maastrichtian faunal 
assemblages^^. This taxonomically standardized and globally 
representative dataset is the same one that was previously used to 
show a geographically uniform pattern of extinction across the K-Pg 
boundary^°. 

Results 

Biogeography above the species level. The ten biogeographic units 
(BUs) we identified have boundaries that are naturally delimited by 
the major patterns of geographic ranges of genera, and reflect sudden 
biogeographic transitions in the biota (Fig. 1). The sampling intensity 
(number of localities per biogeographic region) do not affect BU 
delineations unless poorly sampled regions contain few or no 
endemic taxa, preventing our identification of a biogeographic 
transition. Australia is perhaps the best example of this - our 
analysis groups the single Australian locality with the mostly 
European non-rudist BU (M2). However, the Australian locality 
contains 17 taxa, 15 of which are found on average in five other 
BUs^\ This suggests that, at least given the genus -level data 
available, Australia had cosmopolitan bivalve fauna that extended 
down from European shorelines. In turn, this indicates that either 

i) biogeographic structure above the species level has little to do with 
Mastrichtian Stage oceanic surface currents, and likely reflects 
complex ancient dispersal patterns and plate tectonics, or that 

ii) additional sampling in Australia, for example, and in under- 
sampled localities, may refine and improve the biogeographic 
structure we identified. 

Overall, the genus-level BUs we identified are not necessarily lim- 
ited to localized continental shorelines (Fig. 1). BUs Mj and M2 are 
the best sampled and are distributed mainly in the North American 
Gulf and southern Europe, respectively. BUs M3 and Mg are char- 
acterized by tropical rudist bivalves (Order Hippuritoida) and are 
found in the Eastern and Western Hemispheres, respectively. BU M4 
is composed of several high-latitude clam families, and is distributed 
across South Africa, South America, New Zealand, and the Asian- 
Alaska land bridge. BU M5 is located in both central South America 
and West Africa. BU Me, distributed along the east Asian coastline, is 
a provincial BU with several Asian endemics. BUs M7 and M9 are 

©M^ 

_ ^ OMg 

Figure 1 | The biogeographic structure of bivalves reveals spatial 
organization above the species level. Points correspond to fossil localities, 
and are colored by BU. For visualization, ten-by-ten degree cells were 
colored by BU only if the cell contained fossil localities from a single BU, 
and cells without fossil localities were colored if they were less than 15 
degrees from a locality. Although some cells may have been uninhabitable 
by marine bivalves, they were colored if they met the above criteria. Grids 
are only for visual aid and were not included in any way in subsequent 
analysis. 



provincial and found along North American coastlines, but both 
contain several genera with large geographic ranges. BU Mjo is small 
and tied to the European shoreline; it comprises a single locality with 
many endemics, and has the lowest extinction percentage. However, 
Mio also has an excellent preservational setting, suggesting that its 
difference from M2 may be partially due to sampling effects. The 
full distribution of genera across BUs is available in the online 
supplement^ \ 

BU range as an extinction predictor. The demarcation of BU boun- 
daries has predictive power for per-taxon survivorship. Previous work 
has indicated that geographic range is the best predictor of genus - 
level survivorship in both background and mass extinction^'^^, and 
presumably buffers taxa from extinction, perhaps because it is 
correlated with environmental breadth^^. However, we found that 
BU range is a better predictor than geographic range, or number of 
BUs a taxon is distributed in. A logistic regression with BU range 
alone was the best model to predict K-Pg survivorship (AIG = 
419.74, P = 10"^), compared with a regression with both (AIG = 
421.67, P = 0.78 for geographic range, P = 0.015 for BU range), and 
geographic range alone (AIG = 425.8, P = 10~^). This finding 
suggests that, at least for bivalves in the end-Gretaceous event, the 
number of BUs a lineage was distributed in has better predictive 
power for survivorship than the number of provincial shorelines a 
lineage was distributed in. These results are robust to the inclusion or 
removal (any combination) of inoceramid and rudist bivalves, the 
first whose extinction preceded the K-Pg boundary, and the second 
whose extinction might be tied to a physiological factor [Ref. 25, but 
see Ref. 7]. 

Geographic-range latitudinal selectivity gradient. Gomparing 
extinction percentage across BUs requires an adjustment for the 
correlation of geographic range with extinction probability, 
because differences in distributions of geographic ranges will bias 
observed per-BU extinction percentage. To correct this, we 
estimated expected per-BU extinction percentage given the per-BU 
distribution of geographic ranges^^, and analyzed the difference 
between expected per-BU extinction percentage given geographic 
range and observed per-BU extinction percentage. The entire 
assemblage of bivalves that occurred in each BU was included in 
the calculation of per-BU extinction unless otherwise noted. 

After adjusting the per-BU extinction percentage, we examined 
how extinction percentage varied by BU geography. Whereas Raup 
and Jablonski^° used nine geographic regions to infer that K-Pg 
extinction intensity was globally homogeneous, we used geographic 
regions determined by the modular structure of the data to show that 
BUs had different adjusted extinction percentages (Fig. 2). This 
adjusted extinction percentage did not vary along a paleolongitudinal 
gradient or with distance to either the bolide impact at Ghixculub, 
Mexico or the Deccan flood basalt volcanism of peninsular India 
(P > 0.05). However, we detected a paleolatitudinal gradient of 
adjusted extinction percentage, such that BUs with higher average 
paleolatitudes had higher survival than expected given the geographic 
ranges of the genera that occurred in those BUs (Fig. 2, P = 0.02, 

= 0.495). This result is robust to the inclusion or exclusion of both 
rudist and inoceramid bivalves, but the figure shown has inoceramids 
excluded. The absolute value of the paleolatitude was used to control 
for differences in sampling between hemispheres. 

Discussion 

Our study refines the debate regarding the major causal hypotheses 
for the K-Pg mass extinction. Some researchers have suggested that 
the boUde impact alone triggered a host of secondary effects, for 
example, a brief but intense global thermal pulse^^ and a dust cloud 
that inhibited photosynthesis^^, which together would have led 
to catastrophic extinction cascades. Others have contended that 




SCIENTIFIC REPORT: | 3 : 1790 | DOI: 1 0.1 038/srepOl 790 



2 



o 
o 



I • 

1 



if 

>: 

i 



n ^ ^ ^ \ \ I 

I 0 10 20 30 40 50 60 70 

O 

Average latitude (degrees) 

Figure 2 | Observed survivorship minus expected survivorship, under a 
model that predicts survivorship solely based on geographic range, 
correlates with latitude (Regression with unequal variances, P = 0.02, 
horizontal and vertical standard error shown). These points are not 
independent because BUs contain some of the same genera, but a Mantel 
test between the Jaccard distance matrix and the extinction differences 
between those BUs indicates no correlation between the BU assemblages 
and their adjusted extinction percentages (P > 0.05). Above zero indicates 
survivorship greater than expected given geographic range. BU numbers 
are shown, and rudist BUs are included. Average latitude of each BU was 
determined from the absolute latitudes of the localities in each BU. This 
result is robust to changes of the model parameter, global extinction risk, 
use of median instead of average latitude, and choice of regression analysis 
(equal or unequal variance). Inoceramids are excluded here, but the result 
is not sensitive to their exclusion. If rudists are excluded the P-value 
increases slightly (P = 0.03). 

additional events, such as flood basalt volcanism in India that 
released massive amounts of sulfur and carbon dioxide and resulted 
in severe environmental perturbations, combined with the bolide 
impact to cause the K-Pg mass extinction^^'^^. Our study does not 
reject these hypotheses, but further constrains the search space of the 
proximal causal agents. Our results imply that either the severity of 
the K-Pg kill mechanism declined with distance from the tropics, that 



Mytilus californianus 




Crassostrea gigas 



Ostrea lurida 



Figure 3 | A bipartite occurrence network. Ostrea lurida, Mytilus 
californianus, and Crassostrea gigas each have a second order relationship 
with each other (co-occurrence). Japan and Washington have a single 
second order relationship (shared Crassostrea gigas) . Both Ostrea lurida 
and Mytilus californianus have a single third order relationship with Japan. 
This does not imply that Ostrea lurida, for example, could occur in Japan, 
but more third order relationships than we would expect due to chance 
with Japan is evidence for occurrence potential, or depauperate 
fossilization. 



higher latitude BUs were more resistant to extinction, or both. In 
turn, proposed causal agents and scenarios must either match the 
decline of latitudinal kill mechanism severity with decreased geo- 
graphic range selectivity at the BU level, or provide paleoecological 
evidence for decreased ecosystem extinction risk with latitude. 
Although neoecological evidence suggests that ecosystems at higher 
latitudes have lower extinction risk^°, presumably because of wider 
abiotic tolerances, we cannot be certain that this relationship applies 
through geologic time given differences in latitudinal diversity gra- 
dients and configurations of continents. Analyses of BU-level extinc- 
tion risk in background intervals are needed to calibrate the effect of 
latitude on ecosystem extinction risk throughout the Phanerozoic. 

Though extinction percentage in the K-Pg extinction event was 
geographically uniform^°, our results reveal a selectivity gradient for 
geographic range, such that extinction vulnerability differed depend- 
ing on latitude. Our analysis suggests that taxa must be studied in the 
context of the higher-level biogeographic structure in which they 
reside. In other words, the extinction vulnerability of organisms 
during the K-Pg was biogeographically coupled. Biogeographic 
selectivity, or a difference between BU vulnerability, in the K-Pg mass 
extinction has theoretical ramifications as well. If BUs have differ- 
ential extinction, then invasives from BUs with lower extinction 
might displace BUs with higher extinction as ecospace opens^^ effec- 
tively acting as colonizers for the BU. Given the spatial complexity of 
the K-Pg recovery^, the ability of BUs to displace one another 
through biological invasions is a possibility. Our study underscores 
that the macroecological context of taxa should not be ignored 
because extinction vulnerability is inextricably coupled to biogeo- 
graphic history. 

Methods 

Paleobiology database download. We downloaded the Maastrichtian Raup and 
Jablonski dataset from the Paleobiology Database on March 21, 2012. This dataset 
contains 3,445 occurrences of 329 bivalve genera in 105 assemblages^*^. These 
assemblages are basin level resolution from the Maastrichtian Stage. We chose this 
dataset because it is spatially well sampled and it is taxonomically standardized. We 
could not use all Maastrichtian marine invertebrates (or even bivalves) from the 
Paleobiology Database because the majority of the taxa are from a USGS data dump. 
These data were not used because they are spatially uneven (Gulf of Mexico bias) and 
taxonomically inconsistent with the rest of the database (causing duplicates of many 
genera). Additionally, we opted to use this dataset to make our results more 
comparable to those of Raup and Jablonski^". 

Determining survivors. Survivors were determined from the Paleobiology Database 
standard stratigraphic range intervals We could not use the Sepkoski 
compendium^^ because it lacked range interval data for 68 genera. 

Data availability and software. Data not immediately accessible in the Paleobiology 
Database are available in the public repository Dryad (http://www.datadryad.org). 
This includes the bivalve network, BU assignments, and the geographic ranges of the 
taxa. Shown in the supplement are the distributions of taxa across BUs, geographic 
ranges, and BU ranges. To infer geographic ranges, we used the province- counting 
approach outlined in Jablonski and Raup^^ Their approach used the biogeographic 
provinces in the Atlas of Palaeobiogeography^^. We created shape files for these 
provinces and used the R package sp to infer geographic range size. 

A network approach for biogeography. Here, we introduce bipartite occurrence 
networks, which contain both localities and taxa as nodes (Fig. 3). The links in this 
network (connections between nodes) are occurrences. This network has convenient 
higher- order biological properties. The set of localities a taxon links to is its 
geographic range, while the number of taxa a locality links to is its richness. A pair of 
nodes that are two links away from each other have a "second order relationship." For 
pairs of taxa, the number of second order relationships is the number of 
co-occurrences, while for localities, the number of second order relationships is the 
number of shared taxa (note that a second order relationship cannot exist between a 
locality and a taxon). Third order relationships in this network are between taxon and 
locality, and appear when a taxon is connected to a locality through an intermediary 
taxon and locality (Fig. 3). 

Classical approaches for biogeographic analysis of occurrence data, such as 
ordination or agglomerative clustering of a locality-locality distance matrix, use only 
second order information^^. Specifically, one cannot reconstruct the geographic range 
of a taxon from the distance matrix used for analysis. A network approach allows one 
to integrate geographic ranges, co-occurrence (taxon-taxon), shared taxa (locality- 
locality), and higher order relationships. The higher order relationships can help 
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biogeographic analysis recover from competitive exclusion or taphonomy - for 
example, aragonitic and calcitic shells may not be preserved together in entire 
biogeographic regions^*^"^*^. For smaller spatial scales, such as in a Pacific Northwest 
intertidal ecosystem, goose barnacles may not occur in the same plot as California 
mussels, yet we would like them grouped within the same biogeographic structure 
(intertidal strip). A network approach can therefore be applied to assemblage data 
collected at any spatial or taxonomic resolution, and has the added advantage that no 
dissimilarity measure is required for cluster analysis of the network^^. 

The taxon-locality matrix M, or occurrence matrix, is a representation of the 
occurrence distributions of taxa (any level) across localities (any spatial scale, from 
plots to counties to provinces). Each entry in the matrix is either one (taxon present in 
locality) or zero (absent) 



1 if taxon j is in locality / 
0 otherwise. 



(1) 



This matrix is the basic representation of a bipartite network, where there are two 
types of nodes: taxa and localities (Fig. 3). The links (occurrence relationships) are 
exclusively between localities and taxa, taxa cannot be linked to taxa, and localities 
cannot be linked to localities. 

Taxa will be connected to one another in complex patterns through the localities 
they occupy, making it difficult to interpret the information in the network without 
first isolating the major patterns. To find biogeographic clusters in the network that 
we wish to study, we must identify these patterns. As biologists, we would like to 
identify the boundaries between biogeographic units, which will be encoded as major 
topological features in the network. Moreover, we would like to capture most of the 
information about specific relationships with broad brushstrokes. 

Network community detection is a methodological process that does exactly this: 
identifies the major topological features of networks^''. Though many algorithms have 
been proposed to do this task'*", the map equation is an excellent candidate for 
biogeography because of its accuracy''". 

The map equation approach minimizes as an objective function the theoretical 
limit of a description length of a random walk on the network, where nodes are 
aggregated into community structures. The following story provides some intuition 
about this process. A scientist looks at a random locality. She then randomly chooses a 
taxon found at that locality. Next, she randomly selects a locality that the taxon she 
chose is found in. She now chooses another taxon from the new locality, and repeats 
this process forever. In a network with emergent biogeographic features, she will likely 
spend long bouts of this process within biogeographic units, which we will refer to as 
BUs. Specifically, she can only switch to a new biogeographic unit when she chooses a 
taxon that is not endemic to the biogeographic unit she is currently in. If she would 
like to communicate a list of the localities and taxa she chose, it would save her time to 
communicate a list of the major biogeographic units (which contain taxa and 
localities) she visited, rather than each individual locality and taxon. 

Here we write the basic form of the map equation in the notation of the occurrence 
matrix. However, for more in-depth discussion and derivation, we refer the reader to 
a map equation tutoriaP\ or the original paper^^ 

This bipartite network is unweighted and undirected, all links are equal, and each 
link is symmetric. We refer to the total localities in the network as A, and the total 
number of taxa as T. In undirected networks, the probability that the scientist visits 
any node (locality or taxa) in her infinite random walk is the number of links that 
node has, divided by two times the number of links in the entire network. We refer to 
the number of links multiplied by 2 in the network as L, defined as 



(2) 



We multiply the number of links by two because each link is symmetric, so we need to 
count each link twice. The frequency with which she visits locality i in her random 
walk is therefore 



Pa-- 



while the frequency with which she visits taxon j is 
Pt— j 



(3) 



(4) 



As cartographers, we would like the scientist to convey her random walks as concisely 
as possible. This is an optimization problem - out of all possible partitions, we must 
choose the one that allows her to convey her random walk the most concisely. A 
partition is a division of the network into BUs, such that each node is uniquely 
assigned to a BU. The best partition will be the optimal compression of the major 
patterns in the bipartite network. To evaluate these partitions, we must express the 
frequency with which she switches between BUs. The frequency with which she leaves 
BU m is the number of links that lead from nodes in BU m to nodes outside BU m, 
written and divided by the total number of links 



while the frequency with which she is inside BU m is 



(6) 



where is the total number of links between nodes in BU m, multiplied by 2 because 
each link must be counted twice. These probabilities are sufficient to express the 
extended form of the map equation'' ' 



(7) 



-2E(pr)iog(pr) 
+ E(p"+pr)i«g(pr+pr) 

r A 
Y.^^tHp't+Y.PaHPa 



where L(M) is the amount of information required to convey an infinite random walk 
on partition M. All logarithms listed are in base 2, because information is measured in 
units of bits^^. The fourth term is independent of the partition M, while the first three 
terms change based on the proposed partition. To seek the best partition among 
many, we used an algorithm that is freely available online at (http://www. 
mapequation.org). An applet that illustrates the concept behind the map equation is 
available online (http://www.mapequation.org). 

Partition robustness. The algorithm that minimizes the map equation returned 
largely equivalent partitions across random seeds (7.8-7.92 bits, 1.3-1.33% 
compression). Choice of partition did not change our results. 

Geographic range null model. We use a modification of the approach outlined in 
Payne and Finnegan^^ to test the relative magnitude of geographic range selectivity 
between BUs. What follows is a "fields of bullets" null model, where each genus is 
exposed to a globally homogeneous but province-level kill probability, Pr(die). We 
use this approach because it generates a correlation between survivorship and 
geographic range. In this model, the genera that are in more provinces have more 
chances to survive an indiscriminate culling. Such an assumption generates a 
correlation between geographic range and survivorship. The probability that a genus 
goes extinct is the probability that it dies in every province 

Pr(extinction|?7,)= Pr(die)"', (8) 

where rii is the number of provinces genus / is in. We presume that genus survival is 
binary, surviving if it evades the kill probability once or more in any combination of 
regions 

Pr(survive|?if) = 1 — Pr(extinction|?if). (9) 

The null model is parametrized by Pr(die), therefore the value must be inferred from 
data. We use the least squares estimate of Pr(die) given the observed survivorship data 



Pr(die) = argmin V (X,- - (1 - Pr(die)"'))^ 

Pr(die) ^ 



(10) 



where X, is the survival outcome of a genus, one for survival, zero for extinction. Here, 
argmin is used to indicate a minimization procedure across Pr(die), such that the 
Pr(die) with the least squared error is chosen. Caution should be used when applying 
this least squares estimate to data that do not have many broadly distributed taxa, for 
if too few broadly distributed taxa are in the data the outcome may be overfit to those 
few broadly distributed taxa. However, the Raup and Jablonski dataset has taxa from 
all geographic range levels, so the model is not overfit. The least squares estimate 
Pr(die) for the Raup and Jablonski dataset was 0.805. While changes in this estimate 
could change the expectation of the model, our latitudinal gradient result is robust to 
changes of this probability. The observed survival percentage Yots for a group of taxa 
is therefore 



-E^ 



(11) 



where m is the total number of taxa under consideration. We can also calculate the 
expected survival percentage 7 from the geographic ranges of taxa in the group under 
consideration (for us, a BU). The random variable 7 will be normally distributed 
when the number of taxa is large enough, allowing us to use the standard normal 
distribution function to calculate the probability that the observed extinction 
percentage was generated by the null model - which assumes a province-level kill 
probability that is the same across provinces. The expected value of Y is 



p = — , 



(5) 



iE(l-Pr(die)"'), 



(12) 



(13) 
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while the variance is 

Var(Y)=5]Var(ix,j (14) 

= iE(l-P'-(die)"')(Mdie)"'). (15) 

In the main text, we used this geographic range null model to compare the observed 
BU extinction percentage to the expected BU extinction percentage, or Yobs versus Y. 
The error bars for Figure 2 from the main text were calculated assuming that the 
observed BU extinction percentage was binomially distributed (the variance above 
plus the observed binomial variance) . 
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